First Law of Geography

Everything is related to everything else, but near things are more related than distant things.

  • Waldo Tobler

This to be wary of

  1. Does this actually need to be on a map?

This to be wary of

  1. Does this actually need to be on a map?
  2. Is the structure in my map just another variable lurking in the background?

Recreating Snow's Map

Step 1

Pull in cholera deaths and pump location data.

Notes:

  • Will involve shapefiles, which can be points, paths, lines (vector), or polygons.
  • Based on ESRI's formatting (they develop ArcView GIS)
  • Uses sp and rgdal packages.

DEMO

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/andrewbray/Dropbox/Teaching/2018-spring/math-241/course-materials/data/snow/", layer: "Cholera_Deaths"
## with 250 features
## It has 2 fields
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/andrewbray/Dropbox/Teaching/2018-spring/math-241/course-materials/data/snow/", layer: "Pumps"
## with 8 features
## It has 1 fields

Cholera data

cholera_locs <- coordinates(cholera)
head(cholera_locs)
##      coords.x1 coords.x2
## [1,]  529308.7  181031.4
## [2,]  529312.2  181025.2
## [3,]  529314.4  181020.3
## [4,]  529317.4  181014.3
## [5,]  529320.7  181007.9
## [6,]  529336.7  181006.0
class(cholera_locs)
## [1] "matrix"
cholera_locs <- as.data.frame(cholera_locs)

library(tidyverse)
ggplot(cholera_locs, aes(x = coords.x1, y = coords.x2)) +
  geom_point()

Recreating Snow's Map

Step 2

Pull in map of London. Superimpose cholera and pump data.

Notes:

  • A map involves data originally collected on the sphere.
  • We need both the basemap as well as the projection that will determine how it is displayed on our 2D screen.
  • We'll use ggmap but there are several other packages to do this step.

DEMO

London basemap

library(ggmap)
m <- get_map("John Snow, London, England", zoom = 17, maptype = "roadmap")
ggmap(m)

ggmap(m) +
  geom_point(cholera_locs, aes(x = coords.x1, coords.x2))

This won't work: units down align.

head(cholera_locs, n = 2)
##   coords.x1 coords.x2
## 1  529308.7  181031.4
## 2  529312.2  181025.2
str(m)
##  'ggmap' chr [1:1280, 1:1280] "#FAF5EA" "#FDF9F2" "#FDF9F2" "#FDF9F2" ...
##  - attr(*, "bb")='data.frame':   1 obs. of  4 variables:
##   ..$ ll.lat: num 51.5
##   ..$ ll.lon: num -0.14
##   ..$ ur.lat: num 51.5
##   ..$ ur.lon: num -0.133
##  - attr(*, "source")= chr "google"
##  - attr(*, "maptype")= chr "roadmap"
##  - attr(*, "zoom")= num 17

Projections

A projection transforms coordinates in 3D space to coordinates in 2D space.

library(maps)
map("world", project = "mercator", wrap = TRUE)
## Warning in map("world", project = "mercator", wrap = TRUE): projection
## failed for some data

map("world", projection = "cylequalarea", param = 45, wrap = TRUE)

Projection types

In general, projections can focus on preserving one of the following (at the expense of the other).

  1. Shape/Angle: Mercator, Lambert
  2. Area: Gall-Peters, Albers

Projection information and other info about the geographic data is monitored via a standardized coordinate reference system (CRS) using a library called PROJ.4.

CRS

library(sp)
proj4string(cholera)
## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"

Standard code formats:

  1. EPSG:4326 (WGS84) used by GPS and Google Earth
  2. ESPG:3857 Used for map tiles in Googe Maps, Open Street Maps, etc.
  3. ESPG:27700 (OSGB 1936) commonly used in Britain

CRS

The basemap is actually in EPSG:4326, so let's project our cholera data into those coordinates.

proj4string(cholera)
## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"
cholera_proj <- cholera %>%
  spTransform(CRS("+init=epsg:4326")) %>%
  as.data.frame()
cholera_proj
##     Id Count  coords.x1 coords.x2
## 1    0     3 -0.1363245  51.51291
## 2    0     2 -0.1362775  51.51285
## 3    0     1 -0.1362473  51.51281
## 4    0     1 -0.1362063  51.51275
## 5    0     4 -0.1361612  51.51269
## 6    0     2 -0.1359313  51.51267
## 7    0     2 -0.1365949  51.51285
## 8    0     2 -0.1364394  51.51282
## 9    0     3 -0.1366706  51.51281
## 10   0     2 -0.1366179  51.51292
## 11   0     2 -0.1367315  51.51287
## 12   0     1 -0.1369572  51.51295
## 13   0     3 -0.1368201  51.51271
## 14   0     1 -0.1367729  51.51266
## 15   0     4 -0.1367316  51.51261
## 16   0     1 -0.1370399  51.51273
## 17   0     1 -0.1370924  51.51265
## 18   0     1 -0.1363187  51.51267
## 19   0     4 -0.1362593  51.51260
## 20   0     3 -0.1362059  51.51255
## 21   0     2 -0.1371564  51.51293
## 22   0     1 -0.1371933  51.51308
## 23   0     2 -0.1374392  51.51289
## 24   0     2 -0.1373643  51.51287
## 25   0     2 -0.1372579  51.51290
## 26   0     1 -0.1371470  51.51313
## 27   0     1 -0.1372026  51.51318
## 28   0     3 -0.1372505  51.51324
## 29   0     1 -0.1372814  51.51317
## 30   0     1 -0.1376339  51.51308
## 31   0     1 -0.1377158  51.51315
## 32   0     1 -0.1377106  51.51299
## 33   0     2 -0.1380105  51.51307
## 34   0     2 -0.1381139  51.51303
## 35   0     1 -0.1384685  51.51279
## 36   0     1 -0.1374886  51.51278
## 37   0     1 -0.1380918  51.51250
## 38   0     1 -0.1377217  51.51246
## 39   0     2 -0.1377120  51.51238
## 40   0     8 -0.1375814  51.51245
## 41   0     2 -0.1374305  51.51252
## 42   0     1 -0.1376037  51.51252
## 43   0     1 -0.1368215  51.51232
## 44   0     1 -0.1370186  51.51237
## 45   0     1 -0.1364902  51.51202
## 46   0     1 -0.1364290  51.51196
## 47   0     4 -0.1363780  51.51192
## 48   0     1 -0.1364593  51.51190
## 49   0     1 -0.1365882  51.51200
## 50   0     1 -0.1362121  51.51187
## 51   0     1 -0.1360503  51.51194
## 52   0     4 -0.1359787  51.51198
## 53   0     1 -0.1360447  51.51186
## 54   0     1 -0.1358441  51.51183
## 55   0     1 -0.1357703  51.51185
## 56   0     1 -0.1357216  51.51181
## 57   0     1 -0.1353742  51.51203
## 58   0     2 -0.1355748  51.51214
## 59   0     1 -0.1354466  51.51218
## 60   0     1 -0.1360897  51.51245
## 61   0     1 -0.1359274  51.51226
## 62   0     2 -0.1358134  51.51227
## 63   0     1 -0.1357621  51.51222
## 64   0     1 -0.1357193  51.51217
## 65   0     2 -0.1359251  51.51240
## 66   0     2 -0.1359561  51.51254
## 67   0     1 -0.1358603  51.51256
## 68   0     2 -0.1357804  51.51258
## 69   0     3 -0.1357002  51.51261
## 70   0     1 -0.1354832  51.51268
## 71   0     4 -0.1353908  51.51270
## 72   0    15 -0.1352533  51.51274
## 73   0     3 -0.1351719  51.51276
## 74   0     4 -0.1350995  51.51279
## 75   0     5 -0.1348878  51.51265
## 76   0     2 -0.1347238  51.51251
## 77   0     1 -0.1348178  51.51241
## 78   0     2 -0.1349176  51.51238
## 79   0     1 -0.1349937  51.51235
## 80   0     1 -0.1350932  51.51232
## 81   0     1 -0.1352129  51.51227
## 82   0     1 -0.1353670  51.51222
## 83   0     1 -0.1347528  51.51236
## 84   0     1 -0.1350242  51.51221
## 85   0     1 -0.1349787  51.51214
## 86   0     1 -0.1348168  51.51220
## 87   0     1 -0.1347394  51.51211
## 88   0     1 -0.1348317  51.51198
## 89   0     1 -0.1347715  51.51194
## 90   0     1 -0.1345911  51.51195
## 91   0     1 -0.1345358  51.51190
## 92   0     2 -0.1344959  51.51185
## 93   0     1 -0.1344244  51.51176
## 94   0     1 -0.1347040  51.51185
## 95   0     4 -0.1343347  51.51148
## 96   0     2 -0.1342521  51.51157
## 97   0     1 -0.1341943  51.51152
## 98   0     4 -0.1341110  51.51146
## 99   0     4 -0.1335133  51.51137
## 100  0     1 -0.1335386  51.51154
## 101  0     4 -0.1337880  51.51174
## 102  0     1 -0.1338034  51.51165
## 103  0     1 -0.1338663  51.51170
## 104  0     2 -0.1341597  51.51206
## 105  0     1 -0.1342655  51.51207
## 106  0     2 -0.1343699  51.51216
## 107  0     3 -0.1344271  51.51222
## 108  0     1 -0.1345094  51.51228
## 109  0     4 -0.1345745  51.51234
## 110  0     1 -0.1344772  51.51237
## 111  0     1 -0.1345337  51.51243
## 112  0     7 -0.1337228  51.51226
## 113  0     3 -0.1335162  51.51233
## 114  0     8 -0.1330388  51.51202
## 115  0     1 -0.1329156  51.51169
## 116  0     1 -0.1333610  51.51170
## 117  0     5 -0.1334922  51.51264
## 118  0     8 -0.1327876  51.51255
## 119  0     2 -0.1328993  51.51266
## 120  0     1 -0.1328308  51.51259
## 121  0     1 -0.1329876  51.51273
## 122  0     2 -0.1330345  51.51278
## 123  0     1 -0.1331027  51.51287
## 124  0     2 -0.1331502  51.51292
## 125  0     2 -0.1336385  51.51297
## 126  0     3 -0.1332914  51.51291
## 127  0     1 -0.1335520  51.51302
## 128  0     2 -0.1337383  51.51297
## 129  0     2 -0.1334574  51.51308
## 130  0     3 -0.1341951  51.51272
## 131  0     1 -0.1341558  51.51267
## 132  0     2 -0.1341344  51.51262
## 133  0     1 -0.1340391  51.51254
## 134  0     1 -0.1339959  51.51250
## 135  0     1 -0.1338955  51.51237
## 136  0     1 -0.1342263  51.51276
## 137  0     3 -0.1344430  51.51295
## 138  0     3 -0.1345338  51.51292
## 139  0     3 -0.1346220  51.51289
## 140  0     3 -0.1333934  51.51208
## 141  0     1 -0.1331870  51.51207
## 142  0     2 -0.1332902  51.51204
## 143  0     1 -0.1333941  51.51201
## 144  0     1 -0.1318768  51.51263
## 145  0     1 -0.1316591  51.51272
## 146  0     2 -0.1316900  51.51264
## 147  0     1 -0.1313274  51.51275
## 148  0     1 -0.1323920  51.51303
## 149  0     1 -0.1324359  51.51312
## 150  0     2 -0.1325501  51.51313
## 151  0     2 -0.1324852  51.51301
## 152  0     1 -0.1326663  51.51331
## 153  0     1 -0.1326142  51.51321
## 154  0     1 -0.1330981  51.51319
## 155  0     1 -0.1331764  51.51332
## 156  0     1 -0.1334045  51.51341
## 157  0     1 -0.1333173  51.51309
## 158  0     1 -0.1332787  51.51352
## 159  0     5 -0.1326063  51.51338
## 160  0     1 -0.1325292  51.51325
## 161  0     1 -0.1327582  51.51356
## 162  0     2 -0.1328407  51.51364
## 163  0     2 -0.1328734  51.51369
## 164  0     1 -0.1330518  51.51372
## 165  0     1 -0.1327612  51.51381
## 166  0     1 -0.1325726  51.51387
## 167  0     1 -0.1325537  51.51385
## 168  0     2 -0.1324627  51.51387
## 169  0     2 -0.1324791  51.51389
## 170  0     5 -0.1322153  51.51401
## 171  0     1 -0.1323162  51.51399
## 172  0     1 -0.1322437  51.51396
## 173  0     1 -0.1321192  51.51399
## 174  0     1 -0.1321392  51.51404
## 175  0     4 -0.1320696  51.51405
## 176  0     1 -0.1319574  51.51408
## 177  0     2 -0.1318604  51.51407
## 178  0     1 -0.1317867  51.51410
## 179  0     1 -0.1328681  51.51533
## 180  0     1 -0.1336530  51.51469
## 181  0     1 -0.1337894  51.51464
## 182  0     1 -0.1344164  51.51431
## 183  0     3 -0.1351983  51.51433
## 184  0     1 -0.1349775  51.51440
## 185  0     1 -0.1340471  51.51399
## 186  0     2 -0.1339720  51.51423
## 187  0     1 -0.1332537  51.51396
## 188  0     1 -0.1330836  51.51394
## 189  0     1 -0.1332123  51.51434
## 190  0     1 -0.1340977  51.51388
## 191  0     1 -0.1339554  51.51389
## 192  0     2 -0.1340428  51.51383
## 193  0     1 -0.1338093  51.51371
## 194  0     1 -0.1339705  51.51371
## 195  0     2 -0.1337508  51.51364
## 196  0     3 -0.1338695  51.51360
## 197  0     2 -0.1346203  51.51385
## 198  0     1 -0.1347223  51.51382
## 199  0     1 -0.1346163  51.51403
## 200  0     1 -0.1345116  51.51406
## 201  0     1 -0.1344245  51.51408
## 202  0     1 -0.1346607  51.51410
## 203  0     2 -0.1348155  51.51407
## 204  0     3 -0.1353288  51.51400
## 205  0     3 -0.1353251  51.51376
## 206  0     1 -0.1351935  51.51378
## 207  0     1 -0.1351738  51.51355
## 208  0     3 -0.1350901  51.51364
## 209  0     1 -0.1351064  51.51345
## 210  0     1 -0.1345174  51.51352
## 211  0     1 -0.1343525  51.51357
## 212  0     1 -0.1342768  51.51359
## 213  0     2 -0.1341818  51.51362
## 214  0     2 -0.1342435  51.51352
## 215  0     4 -0.1344024  51.51349
## 216  0     5 -0.1344929  51.51345
## 217  0     2 -0.1345638  51.51344
## 218  0     5 -0.1338792  51.51331
## 219  0     5 -0.1337681  51.51349
## 220  0     3 -0.1339759  51.51329
## 221  0     3 -0.1340732  51.51326
## 222  0     1 -0.1342086  51.51322
## 223  0     5 -0.1342988  51.51318
## 224  0     4 -0.1343863  51.51316
## 225  0     4 -0.1346111  51.51309
## 226  0     1 -0.1349730  51.51297
## 227  0     4 -0.1350693  51.51295
## 228  0     1 -0.1351586  51.51292
## 229  0     3 -0.1352713  51.51289
## 230  0     2 -0.1353472  51.51285
## 231  0     1 -0.1356248  51.51287
## 232  0     2 -0.1350448  51.51334
## 233  0     1 -0.1348969  51.51337
## 234  0     1 -0.1357609  51.51306
## 235  0     2 -0.1358162  51.51311
## 236  0     3 -0.1358666  51.51323
## 237  0     1 -0.1366945  51.51341
## 238  0     1 -0.1357570  51.51326
## 239  0     4 -0.1363895  51.51299
## 240  0     2 -0.1365337  51.51320
## 241  0     2 -0.1366338  51.51313
## 242  0     1 -0.1366662  51.51320
## 243  0     5 -0.1364772  51.51355
## 244  0     3 -0.1363060  51.51424
## 245  0     2 -0.1361014  51.51428
## 246  0     3 -0.1355020  51.51402
## 247  0     2 -0.1354588  51.51420
## 248  0     1 -0.1368681  51.51180
## 249  0     1 -0.1365176  51.51149
## 250  0     1 -0.1361566  51.51135

ggmap(m) +
  geom_point(data = cholera_proj, 
             aes(x = coords.x1, y = coords.x2, size = Count))

Datum

The datum is the point of reference of a GIS system. To fully specify how to put data on a map, you need both the projection and the datum.

Our original cholera data didn't have a datum specified. Some digging will show that it was meant to be in the British standard, so we set that as follows.

proj4string(cholera) <- CRS("+init=epsg:27700")
## Warning in `proj4string<-`(`*tmp*`, value = <S4 object of class structure("CRS", package = "sp")>): A new CRS was assigned to an object with an existing CRS:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
## without reprojecting.
## For reprojection, use function spTransform

Reproject

cholera_proj <- cholera %>%
  spTransform(CRS("+init=epsg:4326")) %>%
  as.data.frame()

ggmap(m) +
  geom_point(data = cholera_proj, 
             aes(x = coords.x1, y = coords.x2, size = Count))

Add in pumps

proj4string(pumps) <- CRS("+init=epsg:27700")
## Warning in `proj4string<-`(`*tmp*`, value = <S4 object of class structure("CRS", package = "sp")>): A new CRS was assigned to an object with an existing CRS:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
## without reprojecting.
## For reprojection, use function spTransform
pumps_proj <- pumps %>%
  spTransform(CRS("+init=epsg:4326")) %>%
  as.data.frame()

Full map code

ggmap(m) +
  geom_point(data = cholera_proj, 
             aes(x = coords.x1, y = coords.x2, size = Count)) +
  geom_point(data = pumps_proj,
             aes(x = coords.x1, y = coords.x2), size = 3, color = "red", pch = 18)

John Snow's Cholera Map

Let's brush this up by using a basemap with a more similar aesthetic to the original.

m <- get_map("John Snow, London, England", zoom = 17, maptype = "toner")

Homework 7

This homework to be done in pairs (but each in your own repo). Get started.